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We study various properties of bosons in two dimensions interacting only via onsite hardcore 
repulsion. In particular, we use the lattice spin-wave approximation to calculate the ground state 
energy, the density, the condensate density and the superfluid density in terms of the chemical 
potential. We also calculate the excitation spectrum, w(k). In addition, we performed high precision 
numerical simulations using the stochastic series expansion algorithm. We find that the spin-wave 
results describe extremely well the numerical results over the whole density range < p < 1. We 
also compare the lattice spin-wave results with continuum results obtained by summing the ladder 
diagrams at low density. We find that for p < 0.1 there is good agreement, and that the difference 
between the two methods vanishes as p 2 for p — > 0. This offers the possibility of obtaining precise 
continuum results by taking the continuum limit of the spin-wave results for all densities. Finaly, 
we studied numerically the finite temperature phase transition for the entire density range and 
compared with low density predictions. 
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I. INTRODUCTION 

The calculation of the thermodynamic properties of 
dilute Bose gases is an old problem that keeps attracting 
renewed attention. Of particular interest is the ground 
state energy per particle, e(p), where p is the density 
of particles. The three dimensional case, with hardcore 
repulsive interaction, was first treated by BogoliubovEI. 
This was later elaborated by many other authorsB using 
a variety of methods. The main result is that the leading 
term in the ground state energy per particle (in three 
dimensions) is e(p) = Airpa(h 2 /2m), where m is the mass 
of the particle and a the scattering length. For hardcore 
bosons, a is simply the diameter of the particles. This 
result was-recently demonstrated rigorously by Lieb and 
Yngvason.0 The first two corrections to this result have 
also been calculated. 

Clearly, this result cannot be true in two dimensions 
since the units would be incorrect. The units would 
work out if one eliminates the dependence on the scat- 
tering length a, i.e. the hardcore boson diameter. How- 
ever, it is not reasonable to expect the energy to be 
independent of the scattering lengtlx The solution to 
this problem was provided by Schicko in 1971 who used 
the fieM theoretic methods developed and applied by 
Beliaevcl to the three dimensional problem. This method, 
which applies in the very dilute limit, consists of as- 
suming that the dominant contributions come from lad- 
der diagrams which are then summed. Schick found 
that the leading contribution to the energy per parti- 



cle is e(p) = 47rp(?i 2 /2m)|ln(pa 2 )|~ 1 with the first cor- 
rection befog proportional to (ln(pa 2 )) -1 . Hincs et al 
re-analyzedcl the integral obtained by Schick and found 
that the first correction term is larger and is equal ta 
— ln|ln(p)|/|ln(p)| where p = pa 2 /n. Lieb and YngvasontJ 
proved rigorously that the leading contribution to the en- 
ergy is indeed the term calculated by Schick, and that 
the leading correctional is between C(|ln(pa 2 )| _1 ) and 
0(|ln(pa 2 )| -1 / 5 ) which is consistent with the correction 
calculated by Hincs et al. 

As mentioned above, these results are for very small 
densities. 

In the meantime interest in hardcore, basons on two 
dimensional lattices has been on the riseOilil and simula- 
tion algorithms have been improving dramatically. 

In this paper we will study in detail hardcore bosons 
on a two dimensional square lattice. Our approach is 
to do very high precision numerical simulaticps—using 
the stochastic series expansion (SSE) algorithmic^] and 
mean field plus spin-wave analyses. The quantities we 
calculate and measure with simulations are: The den- 
sity, p, as a function of the chemical potential, p, the 
condensate, po(p), the energy, e(p), the excitation spec- 
trum and uj{k). In addition, by twisting the boundary 
conditions, the spin- wave calculation also yields the su- 
perfluid density as a function of p (or p). The numerical 
simulations are done at very low temperature (essentially 
T = 0) mostly on 32 x 32 systems with some 48 x 48 re- 
sults. We study these systems from very low densities up 
to half filling, which is equivalent to studying the entire 
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density range [0, 1] due to the particle-hole symmetry. 

Our goal is to determine the range of densities for 
which the spin-wave and continuum calculations are re- 
liable. In addition we are interested in whether and how 
the spin- wave calculation converges to the continuum cal- 
culation (summation of ladder diagrams). 

The paper is organized as follows. In section II we 
present, in some detail, the mean field and spin- wave cal- 
culations of the various physical quantities at T = 0. In 
section III we compare the zero temperature numerical 
results with the results of section II. In section IV we 
compare the spin-wave calculation with the continuum 
ladder diagram results valid for low densities. In section 
V we briefly examine finite temperatures and finally, in 
section VI we present some conclusions and comments. 



II. MEAN FIELD AND SPIN- WAVES 

We will consider a two dimensional system of hardcore 
bosons on an L x L lattice. The hamiltonian is given by 



H 



a\di) 



fj, 



(1) 



where (ij) denotes near neighbors, cii (aj) destroys (cre- 
ates) a boson on site i, and p is the chemical potential. 
The hopping parameter, t, sets the energy scale and will 
be set equal to one in the simulations. In addition, of 
course, the hardcore constraint should be enforced. This 
means that the number of bosons on a given site, i, as 
measured by the number operator, hi = aja^, can only 
be or 1. This is not convenient to express in terms 
of the bosonic operators since the normal commutation 
relations allow any number of bosons on a given site. 

To enforce the hardcore constraint in a simple way, 
we exploit the exact mapping between this system and 
the two dimensional antiferromagnetic Heisenberg model: 
aj <-* ai <-> , and -hi <-> Sf + 1/2. With this map- 
ping the hardcore boson Hubbard hamiltonian becomes 



H 



(ij) 



N, (2) 



where N — L x L is the total number of sites. 
We start with the mean field calculation. 



A. Mean field 

Tiie mean field calculation proceeds in the conventional 
wayli-i The lowest possible energy state corresponds to 
all spins down (no bosons): |0) = JT. | and a total 
magnetization of —N/2. Increasing p, to add bosons cor- 
responds to increasing the externally applied magnetic 
field in the positive z direction. In the mean field so- 
lution, we take all spins parallel and making an angle 



8 with the positive z axis. When the system is empty 
8 = 7r, and the (site independent) azimuthal angle is 
for the moment. Increasing the chemical potential (mag- 
netic field) the spins will turn to reduce the angle with 
the z axis. The mean field state vector is then 



Y[(u + vS+)\0), 



(3) 



where u = sin(0/2) is the amplitude for the spin to be 
down, and v = cos(0/2) the amplitude for an up spin. 
The total z spin is S z = Scos{6) = N b - N/2 where 
S is the largest possible total spin (= N/2) and iVf, is 
the total number of bosons. These equations yield the 
relation between the particle density, p = N^/N, and the 
angle 8: 

sm 2 (0)=4p(l~p). (4) 

Several quantities can now be readily calculated as func- 
tions of p or 8. The expectation value of the hamiltonian 
in the mean field state, \ip), can be easily calculated and 
gives the free energy per site, F (recall that we are using 
the grand canonical ensemble), 



F = -t sin 2 (0) - gcos(tf) - J. 



(5) 



2 w 2 

To determine the dependence of the chemical potential 
on the angle 8 we minimize F, dF/d8 — 0. This yields 

cos(0) = t., (6) 

which, when used with Eq. |^, gives the particle density 
in terms of the chemical potential, 

1 ^ H 
P= 2 + %f 
Therefore, F becomes, 

F 



(J) 



-4tp 2 , 



(8) 



In the ground state, the internal energy (per site), E, 
which is what is measured numerically, is related to the 
free energy via the expression 



E 



F + p( P ), 



(9) 



= -Apt(l-p), 

where we used Eqs. ^ and 0. 

The density of particles in the zero momentum mode 
(the condensate) is given by: 

p = i(V|a t (k = 0)a(k = 0)|^), 



N 

hi 

= P(1-P), 



(10) 
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where the operators 5^(k) (a(k)) are Fourier transforms 
of a\ (a,i). 

These results will be compared with the simulation re- 
sults below. First we calculate the spin- wave corrections. 



The term linear in the operators is required to vanish 
yielding, 



cos(#) = — . 
v ' it 



(15) 



B. Spin- waves 



Taiadude spin- wave corrections, we proceed the usual 
wavLliij'B'EI. We first rotate the z-axis to align it with the 
mean field magnetization direction, i.e. (4> — 0,8). This 
is accomplished by the rotations 



Sf = Sf cos(0) + Sf sin(0), 

QV _ C'V 

D i — J i ) 

Sf = -Sf sin(0) + Sf cos(0), 



(11) 



where S^ is the spin vector in the rotated frame. Equa- 
tion H becomes 



H = -2tJ2[ (Sf Sf cos 2 (0) + 2Sf Sf cos(0)sin(0) 



+SfSf sin 2 (0)) +S?S*) -M^(-Sfsin(0) 



+Sfcos(0) -£iV. 



(12) 



Now, the spins in the rotated frame are expressed in 
terms of bosonic operators b\ and bi. These operators 
do not describe the original hardcore bosons, Eq. In- 
stead they describe the low energy bosonic fluctuations, 
the spin-waves. As such, they obey the usual bosonic 
commutation relations, &t] = § s This transforma- 
tion is accomplished by using 



s? = \{s[ + + s[-)= l -(b\ + b l ), 

S? = ^{S' i + -S'r) = l i {i-h), 



(13) 



Sf = --6^. 



Substituting these expressions in Eq. [12], we assume a di- 
lute gas of spin-waves and ignore cubic and quartic terms 
in the bosonic operators. This yields, 



H 




+ (l + cos 2 (0))(&^+6 l 6t^ 
+(4t sin 2 (6>) + (jt cos(0)) b\bi 

i 

-jv(|(l + cos(0))+t sin 2 (0)) 

+ (-2t cos(0) + |)sin(0)^(bt + b, 



(14) 



This is the same relation between and fi we found in 
the mean field case: Spin-waves do not modify it. 

We are, therefore, left with a quadratic hamiltonian 
which we can diagonalize by first going to Fourier space, 



-ik.r? t 



: lkr 6 k . 



(16) 



This give, 



H = H 



MF 



E(^( & k&k+&l k &-k) 



--BfcC&i&U + fckb-i 



(17) 



where, 



+ -4), 

ft = l(i-(|)V, 

7k = COS^) + cos(fc y ). 



(18) 



Hmf is the constant term in Eq. ful w hich is identical to 
Eq. ||. Eliminating 8 by using Eq. [15| we get the second 
of Eq. ||. In other words, Hmf is just the pure mean field 
result. 

To diagonalize the quadratic spin-wave hamiltonian, 
Eq. [I?], we use the Bogoliubov transformation, 



Ufctt-k, 



(19) 



where a k («k) is a creation (destruction) operator for a 
quasi-particle excitation of momentum k. These opera- 
tors satisfy the bosonic commutation relation, [oik, a k ,] = 
i5k,k', if u\ — v\ — 1. We can ensure this by writing, 

u k = cosh(0 fe ), v k = sinh(^ fc ), (20) 

where <p k is determined by the requirement the Bogoli- 
ubov transformation diagonalize the hamiltonian, Eq. |l7j. 
We, therefore, obtain 



sinh 2 ( 
cosh 2 ( 



1 



A, 



VAT 



(21) 



with A k and B k given by Eq. |l8|. The diagonalized spin- 
wave hamiltonian then becomes, 
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h = (h mf + J2(JaI-bI-a, 



2 -Blfaiak + aKa- 



(22) 



This immediately gives the spin-wave corrected free en- 
ergy per site, Fsw, 



Fsw — -t 



(23) 



where we wrote explicitly the expression for Hmf using 
Eq. |^. Equation |2^ also gives the dispersion of quasi- 
particle excitations, 



w(k) 



2 J A* 



(24) 



It is easy to show, using Eqs. (y_8|2J), that as k — > 0, we 
obtain the linear dispersion characteristic of phonons, 



c(k^0) = 2yi-(£) 2 |k|. 
This gives the sound speed (critical velocity) as 



(25) 



(26) 



Now that we have diagonalized the hamiltonian, we 
can calculate the spin-wave corrections to the quantities 
already obtained by mean field. The number of particles 
kicked out of the condensate (k = mode) due to the 
spin-waves is given by — {ij}\b^bk\ip) ■ Using Eq. 
and the fact that ak\ip) = (no quasi-particle excitations 
in the ground state), we get, 



1 



At 



n ^ = ^7W^M 



(27) 



The spin-wave corrected condensate density is then, 

Po = p (l-p)-±J2 n *- (28) 

Whereas spin-waves did not modify the relation be- 
tween /i and 9, they do introduce a correction to the 
relation between p and p. To find the new expression, 
we simply use p = —dFsw /dp. This gives 



1 



k^o \ V 



Ak — P>k 



- 1 



(29) 



Now that we have determined Fsw an d p{p), we can 
calculate the boson energy per site (internal energy), E, 
using E S w = F S w + p(p)- This gives, 



Esw 



('-<^ 2 )4.i:(\Ar^-A 

k^O \ 



N^At' 



k^O 

A k -B k 



(30) 



C. Superfluid density 

Unlike the above thermodynamic quantities, the super- 
fluid density, p s , requires special treatment of the bound- 
ary conditions. As is well knownt^l, one can relate p s to 
the "spin stiffness" . To accomplish this, we need to com- 
pare the free energy of the system under periodic and 
twisted boundary conditions. This energy difference is 
related to p s . 

In the periodic case, which we treated above, the az- 
imuthal angle, <fi, was taken to be site independent (equal 
to zero). To implement a twist, we take this angle to be 
site dependent, r , and with a constant gradient such 
that the total twist across the system in both the x and 
y directions is 2-7T, 



2n 

T' 



(31) 



where x (y) is a unit vector in the x (y) direction. Con- 
sequently, the rotation of the reference frame to align it 
with the local spin will be site dependent, 

Sf = (Sl x cos(0) + Sf sm(0))cos(&) - Sf sin(^), 

Sf = (sf cos(0) + Sf sin(0)) sin(&) + Sj"cos(&), (32) 
Sf = -Sf sin(0) + Sf cos(0). 

The rotation and diagonalization now proceed exactly 
as before giving, 



li = At cos(#)cos((5</>), 



(33) 



and 



sw - -y • \ 4t cos (^, 

k 

where the T superscript indicates that the boundary con- 
ditions are twisted. A\ and Bj, are given by 



t 



Bl =- 1- 




/' 



At cos(S<p) 

P x2 



7k 



2 V y At cos(<5(/))' 
and 7k has not changed. 



7k, 



(35) 
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We can now calculate AF< 



sw 



r sw 



AF : 



sw 



Al-Bl 



Fsw, 



,m 2 



(36) 



where we made the approximation cos((50) »1- (5(f)) 2 /2. 
The square of the gradient twist can be related to the 
kinetic energy of the supernuidEj by 



(37) 



with v 2 = v 2 +v 2 . Recalling that t — ft 2 /(2m), we obtain 
the superfluid density, 



(38) 



k#0 



It is very interesting to note that at the mean field 
level, the superfluid density and the condensate are iden- 
tical, pa — p s (see Eqs.( p8| , |38| )). As expected, spin-wave 
corrections deplete the condensate, Eq. ^8], but surpris- 
ingly, the su per fluid density is enhanced since the second 
term in Eq. Bq is positive. Furthermore, both our spin- 
wave calculation and numerical simulations (see below) 
show that even at zero temperature (the present case), 
the superfluid density is never truly equal to the den- 
sity nore to the condensate density. In fact, we always 
have p s < p and, of course, p s > po- It is worth noting 
that at the very small densities present in atomic con- 
densates (of the order of 10 -4 ), our results do indeed 
yield p s ~ Pf-p P, which agrees with the experimental 
observations .EJ 

We now compare these results with numerical simula- 
tions. 



III. NUMERICAL SIMULATIONS 

The numerical simulations were done using the SSE 
algorithm.li3ii3 The results we will present here are for 
a 32 x 32 system with periodic boundary conditions and 
inverse temperature (3 = 45. Since we are, mostly, inter- 
ested in the zero temperature properties, we checked the 
finite temperature effects by doing simulations at (3 = 60. 
In addition, we checked the finite size effects by simulat- 
ing 48 x 48 systems. Within our very small error bars, 
the results were identical to the simulations for 32 x 32 
at (3 = 45. We can, therefore, consider these results to 
be valid in the thermodynamic limit at T = 0. 

One of the quantities of prime interest is the superfluid 
density, p s . In the previous section, we calculated p s us- 
ing the relation between the helicity modulus and the 
difference in free energies with periodic and antipcriodic 



boundary conditions.0 Numerically, 
from the winding number Jill 



Ps = 



m (W 2 ) 



is determined 



(39) 



where W is the winding number, D the dimensionality, 
and f3 = 1/ksT. Since the hopping parameter t is related 
to the mass via t — h 2 j (2ma 2 ) , where a is the diameter 
of the hardcore boson (one lattice spacing in our case) , 
Eq. |39| becomes 



pa 



1 (W 2 ) 
2ta 2 D/3 ' 



(40) 



In section V we will consider the finite temperature KT 
transition from superfluid to normal phases. This tem- 
perature from the superfluid density using the universal 
jump conditions, 



2 Ps 



mk B T KT 



(41) 



with m related to t as above. 

The spin-wave results contain momentum summations 
which we did numerically, typically for L = 32, to com- 
pare directly with the numerical simulations. 

In Fig. [l] we show the bosonic energy, E as a function 
of the particle density, p. The circles are the numerical 
results, with error bars smaller than the size of the points. 
The dashed line shows the mean field result, Eq. Q and 
the solid line shows the mean field plus spin-wave cor- 
rection, Eq. [50| The long dashed line will be discussed 
in the next section. We see that agreement between the 
numerical and spin- wave results is excellent for the whole 
range of densities, < p < 0.5, and particle hole symme- 
try extends this agreement all the way to full filling. 

Figure ^ shows the comparison of the condensate frac- 
tion, po, Eqs. ^ and with the numerical results. 
The line convention is as for the previous figure. Once 
again we see that the agreement between spin- waves and 
simulations is remarkable over the entire density range. 
Numerically, the condensate fraction po was calculated 
as the Fourier transformed of the zero momentum equal 
time Greens function. 

In Fig. U we show the dependence of the particle den- 
sity on the chemical potential, Eqs. fj] and Again we 
see that the agreement between the spin- wave calculation 
and the numerical simulation is excellent. 

The superfluid density, p s (Eq. [}8|) is shown with the 
condensate fraction, po (Eq. p9| ), as functions of p in 
Fig. |J. The dashed line is the mean field result for both 
quantities. At the mean field level, po and p s are identi- 
cal which is, perhaps, not so surprising. What is surpris- 
ing is that when we include the spin- waves, po is pushed 
down while p s is enhanced. As we see from the figure, 
agreement with numerical results is extremely good. 

It is clear from Fig. || that for small densities we have 
p s ~ p and po w p. This holds well for p < 0.03 since at 
this density p s and po are no longer equal. 
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FIG. 1. The boson energy, E, vs. the density, p. Circles: 
Simulation results, dashed line: Mean field, solid line: with 
spin- wave correction long dashed: Ladder diagram summa- 
tion of Hines et. alB 
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FIG. 3. The particle density, p, vs. the chemical potential, 
p. Circles: Simulation results, dashed line: Mean field, solid 
line: with spin- wave correction, long dashed: Ladder diagram 
summation of Hines et. alB 



Long dashed: Ladder diagrams 
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FIG. 2. The condensate fraction, po, vs. the density, p. 
Circles: Simulation results, dashed line: Mean field, solid line: 
with spin- wave correction, long dashed: Ladder diagram sum- 
mation of Hines et. ala 



FIG. 4. The dashed line gives both, po and p s vs. p. The 
numerical results are shown in diamonds, p s , and circles, po- 
The solid lines are the corresponding spin-wave results. 
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FIG. 5. The excitation spectrum. The solid line is the 
spin-wave result. It is linear for small |k| with a slope (critical 

velocity) c = 2t^l - (±f . 



As shown in the previous section, the superfluid is sta- 
ble, the low lying excitations are phononic since the dis- 
persion is linear in |k|, Eq. ^5|. In Fig. || we show w(k) 
for p = 0.4. The solid line is from spin-waves while the 
circles are the numerical results. Numerically, w(k) was 
estimated by fitting an exponential to (a(k, r)a^(k, 0)) 
where r is the imaginary time. This quantity is difficult 
to measure£3, especially for large momentum vectors, as 
we see from the noise in the numerical results. It is pos- 
sible to get better numerical results for this quantity by 
using more elaborate methods, such as the maximum en- 
tropy method. 

We see that spin- waves give an excellent description of 
the system. The agreement between the numerical sim- 
ulations and the spin- wave calculation is extremely good 
for all the quantities we considered, except w(k), where 
the difficulty is numerical. Furthermore, this agreement 
extends over the entire density range, < p < 1. 

As seen above, the spin-wave results presented here 
contain momentum sums which we do numerically. We 
can, however, find simple empirical fits which agree with 
the simulations at least as well as the spin-wave calcu- 
lation. The choice of the functional form of the fits was 
guided by two considerations: (a) Particle hole symmetry 
which suggests dependence on p(l — p), and (b) the low 
density diagrammatic results in the continuum (see next 
section). So, we found the following fits for the boson 
energy per site: 



E fit = P (l - P ) I 



47Tfl(l ~ P ) 

IMMi-p))l 



_ \n\\n{A P (l-p))W * 
\HAp(l-p))\J , 



(42) 



where A = 0.537 and B — 3.588 are the fitting parame- 
ters. The condensate density is given by 



fit 

Po 



i 



|ln(0.0335p(l - p))| 
and superfluid density, 

p{ it =p(l-p)(l + 



0.11245 



\Hp(i-p))\ 



(43) 



(44) 



We can obtain an equally good fit for the superlfuid den- 
sity using the form 



Jit _ 



p(l-p)(l + Ap(l-p)), 



(45) 



with A = 0.3356. To give an idea of the goodness of 
these fits, we mention that they agree as well with the 
numerical results as the spin-wave results. It is interest- 
ing to note that, although in Fig. |J p s and po appear 
to have similar functional forms, using Eq. ^5| gives very 
poor results for pq. 



We re-iterate that these fits, Eqs. (|2|j4Jj44j,^5j) are 
purely empirical and are given just as a shortcut to get 
accurate results. 



IV. THE CONTINUUM 

The excellent agreement we found in the previous sec- 
tion between simulations and spin-waves suggests a com- 
parison with the-, ladder diagram summation results in 
the continuum.QiJ The point is that the continuum re- 
sults were obtained using approximations that are valid 
for very low densities, which is where the lattice results 
should approach the continuum. This should give an idea 
of the density where the continuum results start to break 
down. Of course, this breakdown could be due to the 
breakdown of the validity of the ladder diagram calcula- 
tion, or to the increasing importance of lattice effects, or 
both. n 

Using diagramatic methods,J3chickQ showed that the 
condensate density is given by,E3 



1 

Wo 



= -PP 



\f(0,x)\ 2 dx 

X A — \l X 



(46) 



where PP stands for principal part, and p, = 2 pa 2 , po = 
p$a 2 /it. The diameter of the hardcore particles is a. 
In the lattice simulations this corresponds to the lattice 
spacing. /(q',q) is the scattering amplitude in two di- 
mensions for two, particles interacting via the hardcore 
potential. Schicko showed that 



l/(0,q)| 2 = 



16 



[J§(qa)+Y 2 (qa)] 



(47) 



where Jo(qa) and Y^qa) are Bessel functions of zero or- 
der of the first and second kind, respectively. By consid- 
ering the first few terms in the small argument expansion 



7 



of the Bessel functions, Hines et al performed thaintegral 
Eq. Ktq and obtained for the condensate density,B 



-87r 2 lnp 167r 2 ( 7 -ln2)(ln 2 /I-7r 2 ) 



2p /2(ln 2 p + 7r 2 ) 
VSlnV 



/i(ln fi + 7r 2 ) 2 



(48) 



Where 7 is Euler's constant. r-Using this equation with the 
appropriate Green functional, we can obtain the relation 
between the chemical potential and particle density, 

8?rp / ln(-lnp) ( (2 7 + ln2 + 21n7r - 1) 



—lnp \ — lnp 

In 2 (-lnp) 



—lnp 




(47 + 21n2 + 41n?r - 1) 



ln(-lnp) 

(-M 2 



(49) 



Integrating this equation with respect to the density 
yields the energy per particle, 



4np ( ^ m(-lnp) (27 + ln2 + 21n?r - 3/2) 



-lnp \ —lnp 
In 2 (-lnp) 



-lnp 



(-lnp) 2 



(47 + 21n2 + 41n7r - 2) 



In (—lnp) 
(-lnp) 2 



Finally, we can also get the condensate density, 



Pa 



1 



-Hp) 



01 



1 



Un 2 p 



(50) 



(51) 



To compare these equations with the spin-wave and 
simulation results, we should first recall that the kinetic 
energy term used there is the full Laplacian. On the other 
hand, for the lattice bosons we did not include the diago- 
nal term of the Laplacian in the kinetic energy. Further- 
more, we normalized the energy by the number of sites. 
So, the simulation and spin- wave energies should be com- 



pared with p(e p — 4). With this proviso, Eqs. ( [49|[50| 5l[ ) 
are used to calculate the dashed lines in Figs. (|l|^^. 
We see that in all cases the ladder calculation agrees well 
with both the lattice simulations and the spin-wave cal- 
culations for densities < p < 0.05 and sometimes up to 
p = 0.1. It is interesting to note that these densities are 
mucg higher than was thought to be the range of validity 
of the ladder approximation.E 

One can also integrate Eq. 46] numerically thus avoid- 
ing the asymptotic expansion of the Bessel functions. In 
figure |6] we show po versus p from the numerical integra- 
tion of Eq. [l6] (dashed line) and the spin- wave calculation 
(solid line). We see that the exact (numerical) calculation 
of po with the ladder approximation does not improve 




FIG. 6. po vs u. The dashed line is the full numerical 
integration of Eq. Wq, the solid line is the spin-wave result for 
a 1024 x 1024system. 



the comparison with the numerical results. In particu- 
lar, since this approximation is based on the assumption 
of very low density, it does not exhibit the particle hole 
symmetry present in the system. 

To demonstrate how agreement between the ladder 
and spin-wave approximations is approached, we show 
in Fig. the spin- wave energy, Eq. |3^, minus the ladder 
diagram energy, Eq. |5(J versus p. We see that the agree- 
ment is indeed excellent and that the difference tends 
to zero like p 2 . To get a clearer idea of the agree- 
ment, one should also look at the fractional difference, 
AE/Esw = (E sw - Ei adder )/E sw . For example for 
p = 0.022 we find AE/E SW = 1-4 x 10~ 3 , while for 
p = 5x 10~ 5 we find AE/E SW = 1-3 x 10~ 6 . 

While discussing superfluids in the continuum, we take 
the opportunity to make some comments regarding the. 
Onsager-Penrose estimate of the condensate density.EJ 
Often in the literature (see for example reference |2^ 
and references therein) results for the condensate den- 
sity, po, or condensate fraction, po/p, pie compared with 
the estimate by Penrose and OnsagerEil as the standard 
calculation. However, the estimate by Penrose and On- 
sager is incorrect because the trial wavefunction they use 
it too restrictive: It is unity if the particles do not touch 
and vanishes if they do touch. In fact, wavefunctions do 
not vanish so brutally when they approach an infinite 
wall, they go to zero gradually: They start decreasing 
before the wall is reached. It turns out that this effect, 
neglected in the calculation of Ref. 
suits a great dealB 3 ] In three dimensions, 
Ref. ^l], po ~ oil — (47r/3^ 3 p) while using the approach 

DvR O 



changes the re- 
according to 



of BogoliubovEl, one findsE 3 ] po ~ p(l — (8/3)-^(a 3 p)/ir, 
where a is the diameter of the particles. It so happens 
that for the physical values present in experiments, both 
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estimates give po/p) of the order of 10%. However, they 
give very different functional dependence on the density, 
p. Similar differences occur in two dimensions, of course. 



V. FINITE TEMPERATURE 

In addition to the ground state, we studied the finite 
temperature Kosterlitz-Thouless transition to superfluid- 
ity as a function of the particle density. 

To determine the transition temperature, Tkt we ex- 
ploit the universal jump in the superfluid densiy shown 
in Eq. [y]. The results are show in Fig. pi The simula- 
tions were done for L x L systems with L — 16, 32, 64, 96 
and extrapolated to the thermodynamic limit. 

For small, but not too small, particle densities, the 
bosons behave approximately like free particles which 
should give p s cx p. The simulations, Fig. confirm 
this. We get, 



Trh 2 , ,nh 2 

p s7 :— = k B T KT ~ 0.740 2 -—p. 
2 m 2 to 



(52) 



On the other hand, for very low densities, and conse- 
quently very low Tkt-, we should see a crossover to a 
nonlinear regime when the correlation length, £, exceeds 



the mean interparticle diai 
of Fisher and Hohenbergcj 
this limit is given by 

k B T, 



ance, p x l d . Using the results 
, the p dependence of Tkt in 



B-L KT 



AttH 2 p 



2to lnln(l/(pa 2 )) 



(53) 



where a is the scattering length, equal to one lattice spac- 
ing in the case of hardcore bosons. Using equations p3 



and |53|, we get an estimate for the upper bound of the 
density at which the crossover occurs, 



pa < e 



< 10 



-no 



(54) 



A value which is far too small to be tested numerically. 
This gives some limit on the usefulness of the result of 
Fisher and Hohenberg. 



VI. CONCLUSIONS 

By direct comparison with numerical simulations, we 
have shown that spin-waves describe extremely well the 
properties of hardcore bosons, or equivalently spin-1/2 
Heisenberg antiferromagnet in a magnetic field. This 
remarkable agreement extends over the entire density 
range, < p < 1. 

This verification of the accuracy of the spin- wave ap- 
proach in the absence of interactions other than hard- 
core contact, offers a measure of confidence in the ac- 
curacy of such an analysis_even in the presence of near 
and next near interactions .□ We believe a mean field plus 
spin-wave approach should still work well in this case. 
A direct comparison with numerical simulations would, 
however, offer more conclusive assurance. Of course, the 
quality and reliability of the spin-wave calculation de- 
pends strongly on the mean field states that are assumed 
therein lies the challenge. In the case of pure hardcore 
repulsion, the mean field state was easy to guess: Uni- 
form particle density. However, the presence of compet- 
ing interactions, such as near (nn) and next near (nnn) 
neighbors, can complicate things. For example, if the 
nn interaction dominates, it seems reasnoable to use a 
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checkerboard density wave as the mean field state. How- 
ever, there are some other exotic states that might lower 
the energy and need to examined. One such candidate is 
the bond-ordered state.c3 

The agreement between the spin-waves and the lad- 
der diagram approximation for small densities, not only 
bears out the exact leading term of the energy, e(p) — 
4np(h 2 /2m)\\n(pa 2 )\~ 1 , but also agrees very well with 
the logarithmic corrections. For example, Eq. ^0| agrees 
very well with the numerical results, Fig. |l|, up to p rj 
0.12, whereas, keeping only the leading contribution, de- 
viations are observed at p « 0.07. In other words, the 
numerical results are accurate enough to probe the loga- 
rithmic corrections. It would therefore be interesting to 
compare numerical simulations of hard-core bosons in the 
continuum with our spin-wave calculation for large densi- 
ties. We could not find such simulations in the literature. 
One interesting question in that regard is whether the 
continuum limit of the spin-wave equations accurately 
describes bosons in the continuum. In that case, the 
spin-wave approach would offer an alternate simpler way 
to calculate the continuum case. 

Finally, we also determined the KT transition temper- 
ature, Tkt, as a function of particle density. We found 
that the theoretically predicted behaviorci appears to be 
valid at such exceedingly small values as to preclude nu- 
merical and experimental verification. 
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